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FOREWORD 


This  report  was  prepared  by  Physics  International  Company, 
2700  Merced  Street,  San  Leandro,  California  94577,  under 
Contract  No.  F08635-72-C-0229  with  the  Air  Force  Armament 
Laboratory,  Eglin  Air  Force  Base,  Florida.  Mr.  Lovonia  J. 
Theriot  (DLYA)  managed  the  program  for  the  Armament  Laboratory. 
This  effort  was  conducted  during  the  period  from  June  1972  to 
August  1973. 


The  contractor  report  number  assigned  is  PIFR-430. 


This  report  is  divided  into  two  volumes.  Volume  I  presents 
the  generalized  analytical  approach  to  shaped-charge  warhead 
design.  Volume  II  describes  the  modification  and  utilization  of 
a  two-dimensional  finite  difference  continuum  mechanics  code 
utilizing  the  Lagrangian  coordinate  system  to  calculate  the 
complete  jet  formation  parameters  for  any  generalized  axisym- 
metric  shaped  charge.  This  is  Volume  II. 


This  technical  report  has  been  reviewed  and  is  approved. 
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ROBERT  W.  DILLON,  Colonel,  LEAF 
Chief,  Weapon  Systems  Analysis  Division 


ABSTRACT 


This  report  describes  a  technique  to  optimize  the  current 
shaped-charge  design  procedure  as  follows.  Starting  with  the 
desired  target  to  be  defeated,  a  determination  of  the  desired 
penetration  characteristics  of  the  jet  would  be  made*  Exist¬ 
ing  jet  penetration  theory  would  then  be  used  to  estimate  the 
ideal  characteristics  of  the  jet  to  defeat  the  given  target. 

A  shaped  charge  launcher  would  then  be  designed  to  give  these 
ideal  jet  characteristics.  However,  a  suitable  design  pro¬ 
cedure  requires  (1)  a  viable  analytical  or  empirical  design 
approach  to  obtain  a  first  cut  shaped  charge  design,  (2)  a 
better  understanding  than  now  exists  of  the  detailed  mechan¬ 
isms  of  jet  formation,  and  (3)  a  better  understanding  of  the 
phenomenon  of  jet  penetration.  This  report,  which  is  con¬ 
tained  in  two  volumes,  addresses  the  first  two  of  these  re¬ 
quirements.  Volume  I  describes  the  use  of  the  existing  non¬ 
steady  state  theory  of  jet  formation  with  experimental  data 
and  one-dimensional  finite  difference  continuum  mechanics 
calculations  to  obtain  the  linear  collapse  velocity  for  gener¬ 
alized  axisymmetric  shaped  charges.  The  results  of  this  work 
are  then  used  to  obtain  nor-unique  shaped  charge  designs 
which  give  the  required  idealized  jet  parameters.  Volume  II 
describes  the  modification  and  utilization  of  a  two-dimen¬ 
sional  finite  difference  continuum  mechanics  code  utilizing 
the  Lagrangian  coordinate  system  to  calculate  the  complete  jet 
formation  parameters  for  any  generalized  axisymmetric  shaped 
charge.  The  utilization  of  this  code  allows  a  more  detailed 
study  of  such  phenomena  as  jet  stability,  bifurcation  cn  the 
axis,  shear  gradients,  viscosity,  shocks,  incipient  vapor  iza- 
tion,  surface  tension,  and  possible  other  effects.  The  com¬ 
bined  use  of  both  the  engineering  formulations  along  with  the 
sophisticated  two-dimensional  code  calculation  allows  design 
engineers  the  versatility  to  design  the  most  optimum  shaped 
charge  for  their  particular  application. 
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SECTION  I 
rTRODUCTION 


The  design  of  advanced  shaped-charge  warheads  regxxires  a 
detailed  understanding  of  the  mechanisms  of  both  jet  formation 
and  jet  penetration.  Much  of  the  previous  research  on  shaped 
charges  has  been  devoted  to  characterization  testing  of  existing 
designs  rather  than  to  development  of  optimal  designs  with 
maximized  effectiveness  against  a  given  target  or  class  of 
targets.  The  tools  for  developing  such  designs  exist:  advanced 
multi-dimensional,  finite  difference,  continuum  mechanics 
(hydro)  computer  program.  Prior  to  the  effort  funded  under 
this  contract,  an  Eulerian  code  (References  1  and  2)  developed 
by  Systems,  Science,  and  Software  for  Ballistic  Re  jarch  Labora¬ 
tories,  known  as  BRLSC  (Ballistic  Research  Laboratory  Shaped 
Charge) ,  was  the  principal  code  used  to  calculate  the  formation 
of  shaped-charge  jets.  This  code  has  certain  inherent  limita¬ 
tions  : 

(a)  The  HE  burn  cannot  be  calculated  accurately 
without  either  excessive  computer  time  or  cumber¬ 
some  pressure  boundary  condition  patches . 

(b)  The  code  is  limited  to  a  constant  thickness/ 
constant  slope  liner. 

(c)  The  details  of  jet  formation  are  not  well 
described. 
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In  order  to  design  more  effective  shaped  charges  it  is 
~/'essary  to  understand  the  physics  of  jet  formation,  thus  a 
Lagrange  formulation  of  the  problem  is  required.  The  objectives 
of  this  progr  im  were  (a)  to  adapt  the  existing  PJ.SCES  2dl  code 
provide  a  new  Lagrangian  shaped-charge  code  capable  of  cal¬ 
culating  jet  formation  from  variable  thickness,  axially  symme¬ 
tric  shaped  charges,  and  (b)  to  develop  a  generalized  analytical 
theory  to  provide  approximate  geometrical  designs  for  specific- 
application  shaped-charge  warheads. 

The  generalized  analytical  approach  to  shaped-charge  warheac 
design  is  presented  in  Volume  I  of  this  report.  This  volume 
describes  the  modifications  of  the  PISCES  2DL  code  developed  to 
meet  the  objectives  discussed  above  and  presents  the  results 
of  a  computer  simulation  of  the  formation  of  an  axial  jet  from 
the  105  rom  BRL  precision  shaped  charge. 

Section  II  of  this  report  discusses  the  general  calculational 
problems  associated  with  the  numerical  simulation  of  a  shaped- 
charge  jet  and  the  modifications  to  the  standard  PISCES  2DL  code 
that  were  made  to  address  these  problems.  Section  III  presents 
the  calculational  procedure,  and  Section  IV  presents  results, 
conclusions,  and  suggestions  for  further  study. 

The  calculation  was  performed  on  a  CDC  7600  with  170  K 
octal  words.  There  is  no  particular  limitation  on  the  size 
(number  of  zones)  for  a  problem  because  of  the  fcasic  design  of 
PISCES  2DL.  The  modifications  detailed  in  Section  II  were 
straightforward  additions  or  extensions  to  the  structure  of  the 
>de. 


SECTION  II 

TECHNICAL  DISCUSSION 


I4.  GENERAL 

The  numerical  simulation  of  a  shaped -charge  jet  is  challeng¬ 
ing.  A  numerical  study  must  account  for  the  following  phenomena: 
(a)  detonation  of  an  explosive,  (b)  interaction  on  an  explosive- 
metal  interface,  and  (c)  formation  of  a  jet  and  stagnation  region 
due  to  the  collapsing  metal  liner.  The  detonation  logic  employed 
must  result  in  a  detonation  wave  of  the  correct  velocity  and 
magnitude.  The  interaction  of  the  explosive  and  the  metal  liner 
may  result  in  instabilities  that  are  numerical  and/or  physical. 

The  material  properties,  including  the  material  strength  of  the 
metal  liner,  may  be  very  sensitive  in  the  jet  formation  and  in  the 
possible  growth  of  instabilities.  Zonal  resolution  in  the  numeri¬ 
cal  simulation  will  have  a  large  effect  on  the  jet  formation  and 
potential  instabilities.  The  large  distortions  experienced  by 
the  metal  liner  induce  calculational  problems.  Rezoning  tech¬ 
niques  can  be  employed  in  a  Lagrange  code  to  handle  the  distor¬ 
tions  without  sacrificing  a  description  of  the  material  flow. 

The  following  section  discusses  the  various  program  modifications 
made  to  address  the  above  problems. 

2.  ZONING  OP  THE  GRID 

The  initial  zoning  used  to  describe  the  goometry  of  the 
105  mm  BRL  precision  shaped  charge  was  selected  with  the  following 
considerations  in  mind: 
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(a)  A  slideline,  a  computational  slip  surface,  should 
exist  at  the  copper-explosive  interface,  the  copper  liner 
being  the  master  material  and  the  explosive  being  the  slave 
material.  The  motion  of  the  master  material  is  governed  by 
its  internal  stress  gradients  and  by  the  external  pressure 
exerted  on  it  by  the  slave  material.  The  motion  of  the  slave 
material  is  governed  by  its  internal  pressure  gradients  and 
is  constrained  to  slide  along  the  boundary  of  the  master 
material. 

(b)  The  copper  liner  should  contain  at  least  six  columns 
across  its  thickness  and  should  contain  enough  rows  so  that 
no  rectangular  zone  has  a  length-to-width  ratio  greater  than 
two. 

(c)  The  zones  in  the  explosive  should  be  as  square  as 
possible  so  that  the  detonation  phenomena  will  be  described 
correctly. 
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(d)  To  avoid  excessive  tangling  of  the  explosive  slave- 
points  as  they  slide  along  the  copper  liner  (masterline) , 
the  rows  connecting  the  gridpoints  on  the  column  adjacent 
to  the  column  of  slavepoints  (slideline+1)  to  the  slave- 
points  should  intersect  the  mas ter line  at  approximately 
right  angles. 


Several  modifications  had  to  be  made  to  the  standard 
PISCES  2DL  file  to  allow  for  the  above  considerations. 


(a)  The  maximum  number  of  rows  that  a  problem  could  have 
was  increased  from  60  to  145. 

(b)  The  restriction  that  the  column  of  slavepoints  and  the 
slideline+1  column  have  the  same  number  of  rows  was  removed. 

(c)  Non-standard  coordinate-generating  subroutines,  called 
GENLIN,  EQPZON,  and  TRIZON,  were  used  to  generate  the 
coordinates  of  grid  points  in  the  copper  liner  and  the 
explosive. 


The  initial  grid  configuration  is  reproduced  in  Appendix  A. 
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3. 


NEW  DETONATION  LOGIC 


A  new  method  to  simulate  the  ignition  of  the  shaped-charge 
explosive  is  described  in  this  section.  PISCES  1DL  and  2DL  runs 
are  presented  to  show  the  features  of  the  new  method.  In  the 
cne-dimensional  case ,  the  new  method  has  an  advantage  over  the 
standard  method  in  that  the  correct  detonation  velocity  and  sn^pe 
of  the  detonation  front  are  established  in  a  fewer  number  cf  zones 
The  advantage  of  the  new  method  in  the  two-dimensional  case  is 
even  more  pronounced. 

In  PISCES  2DL,  the  standard  high-explosive  burn  logic  is  not 
accurate  for  certain  applications.  Typically,  in  the  standard 
two-dimensional  logic,  the  detonation  velocity  is  correct  but  the 
shape  of  the  wavefront  does  not  agree  well  with  the  shape  predic¬ 
ted  by  theory  or  by  one-dimensional  calculations.  In  particular, 
the  magnitude  of  the  peak  pressure  at  the  detonation  front  is 
typically  only  one-half  to  two- thirds  of  the  Chapman- Jouguet 
pressure  (pCJ)  for  a  reasonable  number  of  zones.  In  contrast, 
the  peak  pressure  at  the  detonation  front  in  the  new  method  is 
equal  to  the  Chapman- Jouguet  pressure,  even  when  the  front  is 
just  passing  over  zones  which  are  adjacent  to  the  zones  being 
detonated. 


a.  Standard  PISCES  1DL  Detonation  Logic .  The  standard 
method  of  describing  an  explosive  detonation  in  PISCES  1DL  is 
presented  in  Appendix  A  of  Reference  3.  To  begin  a  detonation  in 
an  explosive  material,  a  particular  zone  is  specified  by  input 
to  be  "ignited."  Ignition  is  achieved  by  setting  the  zone's  burn 
fraction  initially  to  unity.  This  allows  the  zone's  chemical 
energy  to  be  converted  to  pressure  through  the  zone's  equation  of 
state.  This  pressure  propagates  into  neighboring  zones. 


compressing  them.  Explosive  neighbors  then  will  subsequently 
detonate  if  compression  reaches  the  Chapman-Jouguet  volume.  A 
detonation  in  an  explosive  materia]  will  usually  proceed  for 
15  to  20  zones  in  the  explosive  material  before  the  detonation 
front  is  correctly  established. 

A  problem  illustrating  this  latter  point  is  shown,  ^he 
problem  consists  of  a  5-cm  slab  of  Composition  B  explosive  which  i 
divided  into  100  equally  spaced  zones.  Both  surfaces  of  the  slab 
are  specified  to  be  free  and  the  point  of  ignition  is  specified 
to  be  in  the  leftmost  zone.  The  explosive  is  assumed  to  be 
properly  described  by  a  JWL  equation  of  state  (Reference  3)  with 
the  following  coefficients: 

A  =  5.24  (Mbar)  Rl  =  4.2  u  =  0.34 

B  =  0.0768  (Mbar)  R2  =  1.1 

The  results  of  the  illustrative  problem  are  summarized  in 
Figures  1  through  7.  Figure  1  shows  the  peak  pressure  in  each 
zone.  It  is  seen  that  the  pressure  initially  increases  with 
zone  number  and  then  asymptotically  approaches  a  value  equal  to 
the  Chapman-Jouguet  pressure  of  the  explosive.  This  value  is 
essentially  reached  after  the  explosive  has  burned  through  20 
or  30  zones.  Figures  2  through  7  are  pressure-time  histories  at 
°"'e,  ?1  locations  (zones)  in  the  explosive.  Again,  it  is  ob¬ 
served  that  the  maximum  pressure  obtained  in  a  zone  increases 
with  zone  number. 

b.  Special  PISCES  1DL  Detonation  Logic .  Figures  8  through 
10  show  output  from  a  calculation  similar  to  the  one  described  in 
paragraph  a  (although  it  was  not  carried  as  far  in  time)  using  a 
special  detonation  ignition  logic.  As  can  be  observed  from 
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Zone  number 


Figure  1.  Peak  Pressure  in  Each  Zone  for  Standard 
(  Detonation  in  Composition  B 
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Time  (msec) 


Figure  2 
Position  = 


Stress  Versus  Time  at  Zone  Number  6  (Original 
0.260  cm)  for  Standard  Detonation  in  Composition  B 


Time  (usee) 


Figure  4.  Stress  Versus  Time  at  Zone  lumber  16  (Original 
Position  =  0.760  cm)  for  Standard  Detonation  in  Composition  B 


Time  (nsec) 


Figure  S 
Position  = 


Stress  Versus  Time  at  Zone  Number  21  (Original 
1.01  cm)  for  Standard  Detonation  in  Composition  B 
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Figure  8.  Peak  Pressure  in  Iiach  Zone  for 
Spa- iai  Detoaition  in  Composition  B 
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Figure  9.  Stress  Versus  Time  at  Zone  /'timber  6  (Original 
Position  =  0.26  cm)  for  Special  Detonation  in  Composition  B 


Figure  10,  Stress  Versus  Time  at  Zone  Number  11  (Original 
Position  =  0.S10  cm)  tor  Special  Detonation  in  Composition  R 
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Figures  8,  S,  and  10,  the  new  logic  no  longer  requires  20  or  more 
zones  for  the  detonation  front  to  build  up  to  the  Chapmarv -dengue  t 
pressure. 

The  explanation  for  the  difference  lies  in  the.  maimer  in 
which  the  two  calculations  are  ignited.  In  the  first  case,  the 
ignition  is  accomplished  by  giving  z.  zone  an  initial  burn  fraction. 

war* 

eqbal  to  1.0  and  the  detonation  automatically  propagates  accord¬ 
ing  to  the  burn  logic  described  in  Appendix  A  of  Reference  3. 

In  the  second  case,  the  zone  to  be  ignited  is  assigned  a 
special  equation  of  state  which  determines  the  pressure  in  that 
zone  as  a  pre-programmed  function  of  time.  The  pressure  in  that 
zone  is  forced  artificially  to  rise  from  p  =  0  co  p  «  PCJ  in  a 
time,  x ,  equal  to  the  time  that  it  would  take  for  a  fully  estab¬ 
lished  detonation  wave  to  pass  over  that  zone. 


Except  for  this?  '•  Ifrerence  in  ignition,  the  detonation 
mechanisms  are  identical  in  the  two  sample  calculations.  Each 
uses  the  same  burn  logic  and  the  same  material  description  for 
the  explosive. 

c.  Two-Dimensional  Detonations .  Two  calculations  were 
performed  on  PISCES  2DL,  contrasting  the  standard  two-dimensional 
detonation  logic,  described  in  Appendix  B  (case  A) ,  in  Reference  3 
and  a  special  two-dimensional  detonation  logic  (case  B) .  The 
results  of  these  calculations  are  summarized  in  Appendixes  C 
and  D,  respective ly,  of  Reference  3.  (Figures  C-l  through  C-16 
for  case  A.  are  parallel  in  origin  to  Figures  D-l  through  D-16 
for  case  B.)  The  geometry  for  tne  two  calculations  is  shown  in 
Figure  11. 
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The  burn  logic  used  in  case  B  (special  two-dimensional 
detonation)  is  very  similar  to  the  standard  f  irm  and  involves 
only  slight  modifications  of  the  standard  bur.;  description.  In 
particular : 


(1)  As  in  the  one-dimensional  case  described  earlier,  the 
zone  to  be  ignited  is  assigned  a  special  equation  of  state. 
This  equation  of  state  forces  the  pressure  in  that  zone  from 
0  to  P_T  in  the  characteristic  zone  burn  time,  t- 

(2)  The  burn  fraction  computations  in  the  remainder  of  the 
explosive  zones  do  not  require  a  pre-specif ied  propagation 
velocity,  In  other  words,  the  burn  fraction  calculation 
for  the  special  two-dimensional  calculation  is  the  same  as 
that  described  in  paragraph  b  for  the  one-d‘ mensional 
calculations. 


A  comparison  <■>*=  the  two  calculations  presented  in  Appendixes 
C  and  D  of  Reference  3  shows  significant  differences  in  the  two 
detonation  descriptions.  The  major  distinctions  between  the  two 
solutions  are  that  case  A  appears  to  be  smoother  than  case  B  and 
the  peak  velocities  at  the  detonation  front  of  case  A  are  lov/er 
(by  a  factor  of  1/2  to  2/3)  than  those  of  case  B.  In  addition, 
a  well  defined  detonation  "front"  does  not  exist  in  case  A  but 
does  exist  in  case  B.  Finally,  the  shape  and  magnitude  of 
the  front  in  case  B  appear  to  agree  quite  v/ell  with  theoretical 
predictions  while  the  waveform  in  case  A  is  only  approximately 
correct. 


In  summary,  the  results  of  these  calculations  indicate  that 
the  special  detonation  description  has  a  sharper  and  stronger 
detonation  front  which  agrees  more  closely  with  the  Chapman- 
Jouguet  theory  of  detona cions. 
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4. 


EQUATION  OF  STATE  FOR  COPPER 


Using  data  from  Reference  3,  the  equation  of  state  for 
copper  is  of  the  form 

p  =  (cy  +  Dy^  +  sy  )  x  (l  -  ^  +  r  DE 

where  the  coefficients  have  the  following  values: 

C  =  1.37  (Mbar) 

D  =  1.75  (Mbar) 

S  =  5.6  ((Mbar) 

T  =  1.96 

This  equation  of  state  in  PISCES  is  just  a  form  of  the  polynomial 


equation  of 

state  with 

Al 

=  C  =  1.37  (Mbar) 

A2 

=  (D  ■  ^  C)  =  1.75  -  1.34 

=  0.41  (Mbar) 

A3 

=  (S  ~  -2-  D)  =  5. 6-1. 7  = 

3.9  (Mbar) 

B0 

-  B1  =  T  =  1.96 

To  describe  the  strength  of  the  copper,  a  constant  yield  stress, 
von  Mises  model  is  used.  The  yield  stress  in  uniaxial  tension  is 
0.00219  Mbar  and  the  shear  modulus  is  0.458  Mbar.  The  density  of 

3 

copper  is  8.93  gm/cm  .  A  spall  limit  of  -0.0007  Mbar  was  applied 
to  the  mean  stress  (pressure) . 


SLIDELINE  CALCULATION 


A  more  general  scheme  than  that  presented  in  the  preliminary 
report  has  been  developed  for  the  slideline  calculation.  Calcu- 
lational  problems  can  occur  in  the  motion  of  the  slavepoints  as 
the  explosive  slides  along  the  metal  liner.  A  generalization  of 
the  sliding  logic  allows  for  possible  voiding,  transmission  of  a 
shear  stress  as  well  as  a  normal  stress  across  a  slideline  if 
desired,  and  an  improved  motion  calculation  for  the  slavepoints 
and  the  masterpoints. 

\ 

The  goal  in  the  calculation  of  the  motion  of  gridpoints 
along  a  slideline  is  to  treat  the  calculation  not  in  some  special 
manner  but  in  the  same  way  that  regular  points  are  calculated. 

A  slideline  is  composed  of  a  masterline  and  a  slaveline — two 
separate  adjacent  columns  in  the  Lagrange  grid  (Figure  12) . 
Formerly,  PISCES  2DL  treated  the  motion  of  the  points  on  the 
masterline  (masterpoints)  and  the  points  on  the  slaveline 
(slavepoints)  in  a  special  and  different  way  from  other  points 
in -the  grid.  Howe  rer ,  the  masterline  is  a  geometric  boundary 
for  the  slavepoints  while  the  slavepoints  simply  provide  an 
external  force  to  the  masterline.  A  slavepoint  behaves  precisely 
as  any  other  point  that  exists  on  a  geometric  boundary  and  a 
masterpoint  behaves  exactly  as  a  gridpoint  subject  to  an 
external  stress.  The  logic  used  to  move  a  slavepoint  or  a 
masterpoint  is,  in  theory,  identical  to  that  used  for  any 
boundary  point  under  a  particular  constraint.  Appendix  B 
contains  a  description  of  a  generalized  boundary  motion  scheme 
that  incorporates  masterpoints  and  slavepoints. 


Slave  material 


<ps> 


V  ■  psb 
Figure  12.  Slideline 

As  a  result  of  the  generalized  boundary,  the  weighting  of 
the  external  force  of  the  slave  material  on  the  master  material 
is  done  identically  to  the  way  stress  contributions  around  an 
internal  gridpoint  (Reference  4)  cpntribute  to  its  motion.  Thus- 
rules  regarding  mass  matching  apply  to  the  Lagrange-Lagrange 
interface.  The  mass  weighting  of  stresses  does  not  employ  a  true 
mass  density  times  zonal  volume,  pV,  but  rather  an  areal  density, 
pA,  density  times  zonal  area.  In  planar  symmetry  there  is  no 
difference  between  pA  and  pV  weighting.  In  axial  symmetry  the 
radial  coordinate  of  adjacent  zones  is  taken  to  be  the  same. 

This  is  highly  desirable  across  a  slideline  v'here  a  mismatch  of 
zonal  size  with  non-square  zones  in  conjunction  with  the  radial 
divergence  of  mass  renders  strict  mass  weighting,  pV,  undesirable 
as  well  as  incorrect. 

On  a  slideline  interface  where  different  size  zones  of 
different  materials  may  exist,  it  is  not  important  to  have  exact 
mass  matching.  What  is  more  important  is  to  mass  match  in  a 
normal  direction  to  the  interface.  This  is  identical  to  the  mass 
matching  criterion  of  a  one-dimensional  code. 
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The  code  performs  a  masterpoint  calculation  as  an  interior 
point  calculation  with  the  properly  weighted  surrounding  stress 
contributions.  As  described  in  the  preliminary  report,  a  linear 
weighted  average  of  the  slavepoint  force  contributions  is  utilized 
for  the  force  on  any  particular  masterpoint.  It  should  be  noted, 
mass  weighting  considerations  aside,  that  a  sufficient  number  of 
slavepoints  must  be  defined  along  the  slideline  in  a  tangential 
direction  to  adequately  describe  the  stress  contribution  to  the 
master  material.  This  is  analogous  to  stating  that  sufficient 
points  must  be  used  in  the  piecewise  linear  approximation  of  a 
continuous  function  or,  more  obviously  stated,  that  the  spatial 
resolution  of  the  stress  profile  is  dependent  on  the  number  of 
slave  zones. 

6.  REZONING 

Rezoning  techniques  are  defined  as  methods  whereby  a  dis¬ 
torted  Lagrange  grid  is  mapped  onto  a  new,  smoother  Lagrange 
grid.  This  mapping  may  be  applied  to  grid  coordinates  and  the 
associated  zonal  quantities  of  mass,  energy,  and  stress,  as  well 
as  to  grid  point  velocities.  Rezoning  allows  a  two-dimensional 
Lagrange  calculation  to  be  extended  beyond  a  point  where  the  grid 
has  become  distorted  to  the  extent  that  time  step  inefficiencies 
and  computational  inaccuracies  have  occurred.  Rezoning  is  applied 
at  periodic  intervals  so  that  a  grid  remains  relatively  smooth. 

It  is  important  that  rezoning  not  obscure  (diffuse),  the  solution 
and  yet  allow  a  long  running  time.  The  guiding  precepts  in  the 
use  of  rezoning  techniques  are:  (1)  use  a  rezoner  as  infrequently 
as  possible  and  (2)  when  using  the  rezoner  apply  as  small  a 
perturbation  to  the  solution  as  possible. 


Two  basic  types  of  rezoning  techniques  are  available  at 
the  contractor  facility  for  use  in  conjunction  with  the  PISCES  2DL 
code.  The  first  rezoner  is  known  as  a  manual  snapshot  rezoner 
while  the  second  is  an  automatic  (continuous)  rezoner.  The  snap¬ 
shot  rezoner  is  used  as  follows.  The  user  selects  a  restart  edit 
on  a  restart  tape  to  be  rezoned.  On  special  rezoner  input  cards 
the  coordinates  of  the  new  grid  are  defined.  New  grid  points  can 
be  generated  by  an  equipctential  method,  point  by  point,  or  by 
any  grid  point  generation  technique.  The  rezoner  reads  the 
restart  edit,  generates  the  new  grid,  and  then  maps  the  state  of  the 
old  grid  onto  the  new  grid.  A  restart  edit  of  the  new  grid  and 
rezoned  variables  is  written  on  a  new  restart  tape.  PISCES  2DL 
reads  the  new  restart  tape  and  the  calculation  proceeds  using  the 
new  (rezoned)  grid.  The  snapshot  rezoner  is  integral  to  the 
PISCES  2DL  program,  and  the  rezone  and  restart  of  a  problem  can  be 
run  as  a  single  job. 

For  a  two-dimensional  problem  chat  frequently  gets  into 
calculational  difficulties  due  to  a  grid  distortion,  the  repeated 
application  of  the  snapshot  rezoner  could  become  quite  tedious 
and  time  consuming.  The  automatic  rezoner,  however,  does  not 
require  that  a  2DL  problem  be  stopped  and  restarted  at  rezone 
time.  The  automatic  rezoner  performs  its  rezoning  functions 
simultaneously  as  needed  with  the  normal  calculational  sweep 
through  the  grid.  The  user  specifies  (a)  what  region  he  wishes 
to  be  rezoned,  (b)  what  calculational  cycle  to  start  rezoning, 

(c)  what  cycle  to  stop  rezoning,  (d)  the  frequency  of  cycles  to  be 
razoned,  and  (e)  the  rezoning  options  desired.  Following  the 
general  precepts  of  rezoning,  the  frequency  of  rezone  cycles  is 
kept  to  a  minimum  and  the  re zone  options  specified  perform  a 
small  perturbation  to  the  solution.  It  should  be  noted  that  the 
magnitude  of  the  perturbation  to  the  solution  as  a  result  of  the 
application  of  either  rezoner  is  largely  empirical  and  can  only 
be  judged  from  previous  tests  and  from  the  quality  of  the  solution. 
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The  use  of  one  of  the  two  rezoners  does  not  preclude  the  use 
of  the  other.  In  fact,  the  rezoners  may  be  used  together  and  do 
not  interfere  with  each  other.  A  more  detailed  description  of 
the  snapshot  rezoner  is  contained  in  Reference  5.  Appendix  C 
contains  a  description  of  the  automatic  rezoner. 

7.  ANTI-HOURGLASSING  OPTIONS 


*  Two-dimensional  Lagrangian  codes  that  use  quadrilateral 
zones  are  susceptible  to  a  grid  instability  known  as  hourglassing. 
The  anomalous  hourglass  shaped  grid  distortion  results  in  in¬ 
correct  computed  displacements  and  eventually  more  serious  prob¬ 
lems. 


PISCES  2DL  has  the  option  of  employing  a  triangle  viscosity 
to  damp  the  anomalous  motion  of  the  grid  points.  The  triangle 
logic  utilizes  a  Navier-Stokes  type  viscosity  defined  for  tri¬ 
angular  zones  surrounding  each  grid  point.  The  viscosity 
coefficient  varies  for  particular  materials  and  with  the  amount 
of  numerical  damping  to  be  imposed  on  the  solution.  Reference  7 
gives  a  further  explanation  of  this  technique. 

Another  approach  to  eliminating  hourglass  instability,  known 
as  hourglass  subtraction,  has  been  developed  by  the  contractor. 
The  idea  is  to  subtract  the  hourglass  motion  from  the  grid  at 
each  cycle.  The  hourglass  component  of  the  velocity  fieli.  of  a 
mesh  of  quadrilaterals  is  defined,  the  hourglass  component  having 
zero  net  momerefcwn.  Subtracting  «fchis  hourglass  component  from  the 
velocity  field  at  each  cycle  maintains  a  mesh  of  parallelograms 
and  results  in  a  continuous  velocity  smoothing  which  is  similar 
to  the  velocity  smoothing  option  of  the  continuous  rezoner 
(Reference  7)  . 
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IMP  OPTION 
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Occasionally  in  a  calculation,  a  zone  far  from  the  region 
of  interest  may  distort  to  the  point  where  it  controls  the  time 
step  of  the  entire  problem.  An  example  of  this  occurs  on  the 
free  surface  of  the  unconfined  explosive  in  the  shaped  rharge 
under  study.  The  explosive  blows  off  radially  and  its  turbulent 
nature  causes  the  Lagrange  grid  at  the  free  surface  in  some  aireas 
to  get  into  time-step  difficulty  or  to  tangle  into  a  zero  or 
negative  volume.  Rather  than  simply  applying  the  rezonar  to  this 
region,  it  is  much  easier  to  simply  ignore  these  zones.  This  is 
possible  when  these  zones' liave  released  all  their  chemical  poten¬ 
tial  energy  and  in  fact  can  no  longer  influence  the  motion  of  the 
liner.  Removal  of  these  zones  from  the  problem  requires  that 
they  be  isolated  from  the  rest  of  the  grid.  A  logic  called  IMP, 
Isolated  Mass  Point,  is  used.  The  procedure  consists  of  replacing 
the  zone  with  a  free  boundary  condition  and  not  allowing  the  zone 
to  be  calculated  any  more,  but  rather  maintaining  the  values  it 
possesses  at  the  time  it  is  removed  from  the  grid. 

9 .  TRACER  POINTS 

Implicit  in  the  idea  of  Lagrange  rezoning  is  the  fact  that 
original  zones  and  their  associated  values  will  lose  their 
identity  through  rezoning.  As  rezoning  is  performed,  it  becomes 
more  and  more  difficult  to  tell  where  mass  originated.  This  is  a 
common  problem  for  Euler ian  codes  where  mass  units  never  exist. 
Eulerian  codes  often  use  Lagrange  tracer  points  to  define 
material  interfaces  and  allow  general  mapping  of  the  fluid  motion* 
These  tracer  points  are  moved  through  the  grid  according  to  the 
average  velocity  of  the.  surrounding  material. 
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A  scheme  has  been  devised  whereby  tracer  points  could  be 
used  in  a  Lagrangian  code  in  conjunction  with  a  rezoner.  It 
should  be  obvious  that  tracer  points  are  not  necessary  in  a 
Lagrange  code  if  no  rezoning  is  performed.  Tracer  points  in  2DL 
are  different  than  Eulerian  tracer  points  in  that  their  motion  is 
described  completely  by  the  motion  of  the  Lagrangian  cell  in  which 
they  reside.  A  Lagrangian  tracer  point  stays  with  a  Lagrangian 
zone  and  moves  as  the  zone  moves,  maintaining  its  same  relative 
position  within  the  zone.  At  rezone  time  the  tracer  point  may 
find  itself  in  a  new  relative  position  inside  the  zone  or  in  fact 
may  now  reside  within  a  different  zone.  The  tracer  point  only 
has  to  be  updated  at  rezone  times  rather  than  every  cycle  as  in 
an  Eulerian  code.  The  tracer  point  positions  and  velocities  are 
derived  / from  two  numbers  defining  the  relative  position  inside 
the  zone  and  the  zone's  coordinates  and  velocities.  It  is  only 
necessity  to  derive  the  tracer  positions  and  velocities  at  output 
times  for  printout  and  plotting.  This  scheme,  although  designed, 
was  not  implemented  into  the  PISCES  2DL  code  for  these  calcula¬ 
tions  . 
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SECTION  III 

CAI.CULAT10NAL  PROCEDURE 


The  shaped- charge  calculation  proceeded  as  follows : 


(3)  The  grid  was  generated,  plotted,  and  checked. 


(2)  Calculation  was  begun  at  time  zero  with  initiation  of 
the  explosive  and  run  to  a  point  just  before  the  detonation  wave 
hit  the  copper  liner  5.6  ysec) .  This  duplicates  exactly  the 
previous  solution  (Reference  3) . 


(3)  At  this  time,  on  the  basis  of  the  prelim' nary  calcula¬ 
tion,  it  was  decided  to  employ  the  automatic  rezoner  in  the  apex 
area  of  the  copper  liner.  The  following  automatic  rezone 
options  were  employed: 


Three  rezone  areas  were  defined. 


Region  1 — columns  2  to  6,  rows  120  to  146 
Region  2 — column  1,  rows  120  to  146 

Region  3 — columns  2  to  7,  rows  120  to  146.  (Column  1 
is  the  inside  surface  of  the  liner  while  column  7  is 
the  liner  surface  in  contact  with  the  explosive.) 


The  options  for  each  rezone  region  were  as  follows  (see 
Appendix  C) : 


Region 

1 

NSTART 

~ 

199 

NEND 

=r. 

10000 

NFREQ 

- 

10 

N50PT 

14200 

Region 

2 

NSTART 

— 

198 

NEND 

=: 

10000 

NFREQ 

= 

10 

N50FT 

= 

12.:00 

Region 

0 

NSTART 

= 

198 

NEND 

— 

10000 

NFREQ 

= 

10 

N50PT 

12250 

The  automatic  rezoning  schedule  is  as  follows:  At  cycle 
198  (—6.8  fisec)  ,  column  1/  the  inside  surface  of  the  liner  has 
a  coordinate  rezone  of  NXY  =  2  (Appendix  C) .  This  helps  keep 
the  points  on  this  column  from  crossing  over  each  other  as  they 
form  the  jet.  No  velocity  rezoning  or  momentum  conservation  is 
performed.  The  velocity  distribution  along  the  inner  surface  is 
quite  important  so  that,  rather  than  the  true  phenomenon  being 
obscured,  no  velocity  rezoring  is  performed. 

At  cycle  198,  region  3  has  a  coordinate  rezone  of  NXY  =  2. 
This  will  help  keep  zones  from  tangling  by  trying  to  maintain 
parallelogram  zones.  A  velocity  rezone  of  NVEL  =  5  is  performed 
to  smooth  slightly  the  velocity  distribution  in  the  interior  of 
the  liner  and  the  explosive  interface.  This  small  amount  of 
smooching  is  done  in  order  to  damp  out  the  ripples  previously 
experienced  at  the  explosive  liner  interface.- 


At  cycle  199,  the  next  cycle,  region  1  has  a  coordinate 
re2one  of  NXY  =  4.  This  allows  the  columns  within  the  liner  to 
remain  equidistant  in  order  tc  keep  adequate  resolution  near  the 
expanding  jet.  No  velocity  rezoning  is  performed. 

The  sequence  of  automatic  re zoning  will  then  occur  every 
ten  cycles  thereafter  until  the  user  respecifies  the  options. 

In  general,  a  coordinate  re zone  of  NXY  =  4  should  be  used 
only  after  a  coordinate  rezone  of  NXY  =  2,  because  the  centroid 
method  (NXY  =  3# 4)  will  not  work  properly  on  a  tangled  zone. 

The  operational  rule  is  to  use  NXY  =  1  or  2  to  untangle  the  grid 
and  to  use  NXY  =  3  or  4  to  more  evenly  space  zones .  The  cen¬ 
troid  method  (NXY  =3,  4)  was  used  in  the  case  of  expanding 
zones  to  maintain  good  resolution.  In  the  case  of  compressing 
zones,  NXY  =  3  or  4  will  keep  a  zone  from  crushing  with  a 
resultant  small  time  step. 

Note  that  in  the  three  rezone  regions  specified,  care  was 
taken  in  the  selection  of  rezone  parameters  to  avoid  overly  con~ 
straining  the  solution.  Coordinate  rezonings  vised  the  lesser 
constraint  of  NXY  =  2  or  4  versus  NXY  =  1  or  3.  Velocity 
rezoning  was  only  performed  in  the  interior  of  the  copper  liner 
and  on  the  back  surface.  The  value  of  NVEL  -  5  is  an  inter¬ 
mediate  value  and  does  not  force  a  large  constraint  on  the  solu¬ 
tion  to  the  problem.  Also,  it  is  important  to  note  that  the 
rezoning  does  not  destroy  the  various  gradients  that  exist  in 
the  problem.  Radial  or  axial  velocities  are  not  forced  to 
behave  in  some  monotonic  manner.  This  .is  highly  important  if 
the  actual  jet  formation  process  is  to  be  characterized. 


(4)  Snapshot  rezoning  of  slavepoints  was  performed  to  allow 
for  better  resolution  of  the  stress  profile  along  the  2iner, 

The  slavepoints  wore  moving  in  a  more  reasonable  manner  due  to 
the  new  motion  logic;  however/  due  to  the  relatively  coarse 
explosive  zoning,  tnese  snapshot  rezones  were  necessary. 

(5)  The  calculation  was  carried  to  10  microseconds.  At 
this  point  a  large  snapshot  rezoner  would  have  to  be  performed 
in  the  liner.  This  would  move  the  columns  toward  the  front  of 
the  jet  co  al-ow  for  more  resolution  in  the  jetting  area  and 
would  result  in  the  stagnation  region  having  relatively  coarser 
zoning.  This  is  allowable  since  significant  gradients  do  not 
exist  in  this  region. 
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SECTION  IV 
RESULTS 


The  shaped-charge  jet  calculation  is  presented  in  coordi¬ 
nate  and  velocity  vector  plots  in  Appendix  A.  The  results  from 
the  preliminary  calculations  are  also  presented  in  Appendix  A 
for  comparison.  A  comparison  of  the  previous  and  current  calcu¬ 
lation  at  about  8  ysec  shows  that  the  explosive-metal  interface 
is  experiencing  instabilities  near  the  liner  apex  in  both  calcu¬ 
lations.  In  the  first  calculation,  these  instabilities  con-  ■ 

tinued  to  grow  and  it  was  necessary  to  use  the  snapshot  rezones 
to  eliminate  them  from  the  problem.  Since  the  cause  of  these 
instabilities  was  uncertain,  it  was  decided  to  let  the  auto¬ 
matic  rezoner  try  to  damp  them  out.  As  seen,  the  automatic  . 

rezoning  did  not  get  rid  of  the  instabilities  completely  but  ; 

did  retard  their  growth.  ! 

I 

There  are  twc  possible  sources  of  these  instabilities.  The  ) 

f 

first  may  be  a  characteristic  of  the  slavepo.int  position,  and  the  I 

second- -^ay--be~a™Tay4^r --ins-feab-ii-irty-; — The~posTt±oft'  Of  the  slavepoint 
with  respect  to  the  master  material  would  give  characteristic 
wavelengths  of  the  instabilities  dependent  UjDon  th^ slavepoint 
position.  If  this  is  the  cause  of  the  instabilities,  then  finer 
zoning  in  the  explosive  along  the  liner  would  remedy  it.  It 
should  also  be  noted  that,  prior  to  the  use  of  a  linear  interpo¬ 
lation  scheme  for  the  slave  zones,  very  large  amplitude  insta¬ 
bilities  were  produced  on  the  interface. 
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The  rippling  on  the  explosive-metal  interface  may,  in  fact, 
be  due  to  a  classic  Taylor-like  instability.  The  calculation  at 
the  interface  involves  the  relatively  constant  pushing  of  a, heavy 
material,  copper,  by  a  light  material,  Composition  B,  which  is  a 
condition  ripe  for  a  Taylor  instability.  Thp  ripples,  if  they 
are  a  Taylor  instability,  should  be  worse  in  areas  where  there 
is  a  shorter  radius  of  curvature.  The  observation  is  made  that 
the  ripples  occur  near  the  apex  of  the  cone  and  are  not  par¬ 
ticularly  evidenced  in  the  linear  part  of  the  liner. 

If  the  ripples  on  the  back  surface  of  the  liner  are  Taylor 
instabilities,  the  amplitude  of  the  perturbation  will  decay 
exponentially  with  distance.  For  the  given  width  of  the  liner 
the  amplitude  at  the  back  surface  will  have  decayed  to  about 
5  percent  of  its  value  when  it  reaches  the  inside  surface.  How¬ 
ever,  the  automatic  rezoner  through  successive  rezones  will 
couple  the  motions  on  the  back  and  front  of  the  liner  and  there¬ 
fore  help  transmit  any  perturbation  that  occurs.  It  is  there¬ 
fore  desirable  to  keep  the  ripples  to  a  minimum  in  the  numerical 
simulation.  It  should  be  noted  that  the  Taylor-like  instabili¬ 
ties  may  also  occur  physically  and  be  highly  sensitive  to 
machining  tolerances. 

Independent  of  the  cause  of  these  instabilities  there  are 
means  that  could  be  employed  to  keep  them  from  growing.  The 
first  is  to  use  the  triangle  viscosity  described  in  Section  III, 
and  the  second  is  to  increase  the  copper  yield  strength.  Since 
copper  is  strain-rate  dependent,  it  is  probable  that  the 
2.19-kbar  value  is  too  low  for  the  high  strain  rates  of  this 
problem  (Reference  8) . 
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At  8.7  ijsoc  on  the  current  calculation,  a  perturbance  is 
noted  on  the  liner  that  was  not  associated  with  the  previous 
instabilities.  It  is  believed  that  this  perturbance  is  the 
result  of  an  adjustment  in  the  automatic  rezoner  near  this 
time.  Further  calculations  would  be  needed  to  clarify  this 
assumption. 

The  sequence  of  plots  illustrates  the  nature  of  jet  break¬ 
out.  At  early  times  it  is  noted  that  the  jet  does  not  have  a 
monotonic  character.  The  first  few  zones  just  off  the  axis 
initially  lag  the  rest  of  the  liner  motion.  This  appears  to  be 
the  result  of  the  instabilities  affecting  the  inside  of  the 
liner.  As  the  solution  progresses,  the  jet  can  be  seen  to  form 
in  a  smooth  manner.  The  automatic  rezoner  helps  allow  for  the 
smooth  nature  of  the  solution.  The  application  of  the  automatic 
rezoner  has  shown  itself  to  be  a  highly  useful  tool.  Further 
study  and  experimentation  are  necessary  to  develop  its  potential 
fully.  Additional  algorithms  for  coordinate  and  velocity 
rezoning  are  being  developed  to  increase  its  capabilities. 

Unfortunately,  the  scope  of  this  contract  did  not  permit 
running  the  final  calculation  long  enough  for  complete  compari¬ 
son  with  experimental  data  ana  further  comparison  with  the  pre¬ 
liminary  calculation  was  not  possible.  However,  one  interesting 
result  of  the  preliminary  calculation  was  the  forrr  tion  and 
growth  of  the  jet  tip  instability  followed  by  a  stable  jet. 

This  same  jet  tip  instability  has  been  observed  experimentally 

* 

with  X-ray  photography.  The  jet  velocity  increases  from  about 
0.6  cm/ysec  to  0.78  cm/ysec  from  tne  tip  through  the  instability 
at  a  time  of  13  ysec. 

*Personal  communication  between  Dr.  C.  Godfrey  and  L.  Behrmann 
with  R.  Vitaii  of  BRL. 
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SECTION  V 
CONCLUSIONS 


The  calculation  of  a  shaped- charge  jet  with  a  finite  dif¬ 
ference,  Lagrangian  coordinate  continuum  mechanics  code  has  been 
demonstrated.  Although  the  final  calculation  was  not  carried 
out  as  far  as  the  preliminary  calculation,  a  sufficient  jet 
length  was  obtained  to  arrive  at  the  above  conclusion.  The  fact 
that  che  jet  tip  showed  a  mass  buildup  and  instability  similar 
to  experiments  could  result  only  if  the  jet  radial  velocities 
were  not  arbitrarily  fixed  to  allow  the  problem  to  run. 

The  code  modifications  that  were  made  assured  that  the 
physics  of  the  problem  was  correctly  calculated.  Two  of  those 
modifications  are  worth  reviewing.  The  first  is  the  automatic 
rezoner.  This  rezoner  was  designed  with  considerable  flexi¬ 
bility  so  that  the  engineer  could  adjust  both  the  magnitude  and 
frequency  of  the  rezone  within  a  specified  space.  With  this  and 
the  snapshot  rezoner,  Lagrange  calculations  will  be  more  versa¬ 
tile  where  large  grid  distortions  occur.  The  second  modification 
of  importance  is  the  availability  of  multiple  column  slidelines 
with  the  option  of  void  open  and  closure  on  each  slitleline.  This 
modification  will  allow,  for  example,  the  study  of  possible 
bifurcation  of  the  jet.  Further  modification  of  the  automatic 
rezoner  to  include  tracer  prints  would  be  desirable  but  was 
beyond  the  scope  of  this  contract.  When  rezoning  is  performed, 
correlation  with  the  jet  parameter  theory  developed  in  Volume  I 
of  this  report  would  be  tedious  unless  tracer  points  are  used. 
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In  summary,  this  calculation  addressed  the  solution  of  the 
complete  shaped-charge  jet  formation  problem.  A  full  explosive 
detonation  was  calculated  for  the  loading  on  the  metal  liner. 

The  code  is  not  restricted  by  the  shaped-charge  geometry  as  long 
as  an  axis  of  symmetry  is  present,  i.e.,  variable  thickness 
liner,  curved  liner,  multi-material  liner,  and  explosive  case 
shaping  are  all  acceptable  geometries .  The  code  can  also  use  a 
forcing  function  to  simulate  the  explosive  loading  on  the  copper 
liner.  The  forcing  function  could  be  an  external  pressure  or 
velocity  constraint.  The  velocity  condition  would  not  allow  for 
any  Taylor-like  instabilities  to  appear  or.  the  back  surface  cf 
the  liner.  f 


REFERENCES 


M.  L.  Gittings,  R.  I.  Sedgwick,  and  J.  M«  Walsh,  Numerical 
Analysis  of  Jet  Formation  from  Lined  Shaped  Charges ,  Con- 
tract  Report  No.  51,  3SR-641,  Ballistic  Research  Laboratory 
August  1971. 

L.  I.  Walsh,  D.  E.  Wilkins,  and  R.  T.  Sedgwick,  User's 
Manual  for  an  Eulerian  Shaped  Charge  Computer  Program,  Air 
Force  Armament  Laboratory  Technical  Report  AFATL-TR-73-24 , 
February  1973. 

R.  Hofmann,  N.  Birnbaum,  and  S.  Jardin,  Preliminary 
Calculation  of  Jet  Formation  from  the  105-mm  BKL  Precision 
Shaped  Charge,  TCAM  Technical  Memo  72-21,  Physics  Interna¬ 
tional  Company,  November  1972. 

PISCES  2DL,  Manual  A,  General  Description  and  Finite- 
Difference  Equations ,  Revision  1,  Physics  International 
Company,  San  Leandro,  California,  19/2. 

S.  L.  Hancock,  A  Lagrange  Rezoner,  TCAM  Technical  Memo  73-1 
Physics  International  Company,  San  Leandro,  California, 
February  1973. 

E.  D.  Giroux,  HEMP  User's  Manual,  UCRL-51079,  Lawrence 
Livermore  Laboratory,  June  1971. 

S.  L.  Hancock,  An  Hourglass  Subtraction  Procedure,  TCAM 
Technical  Memo  73-6,  Physics  International  Company,  San 
Leandro,  California,  May  1973. 

E.  L.  Lee,  H.  C.  Hornig,  and  J.  W.  Kury,  Adiabatic 
Expansion  of  High  Explosive  Detonation  Products,  UCRL-50422 
Lawrence  Radiation  Laboratory,  Livermore,  California, 

May  2,  1968. 


39 

(The  reverse  of  this  page  is  blank) 


APPENDIX  A 

COORDINATE  AND  VELOCITY 
VECTOR  PLOTS 


i 

f 


I 


Figures  A-l  through  A-12  are  coordinate  and  velocity  vector 
plots  of  the  preliminary  shaped-charge  jet  calculations.  Figures 
A- 1.3  through  A-33  represent  the  current  calculations. 
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Figure  A- 3..  Snapshot  Velocity  Vector  Plot 
at  Time  =  10.3  ysec 
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Figure  A-4.  Snapshot  Velocity  Vector  Plot 
at  Time  =11.0  ysec 
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Figure  A- 5.  Snapshot  Velocity  Vector  Plo 
at  Time  «  12.3  usee 


51 

(The  reverse  of  this  page  is  blank) 


Figure  A- 6.  Snapshot  Velocity  Vector  Plot 
at  Time  =13.0  nsec 
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Figure  A-7 .  Snapshot  Velocity  Vector  Plot 
at  Time  =8.01  usee 
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Figure  A-ll.  Snapshot  Velocity  Vector  Plot 
at  Time  =  12.3  user 
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Figure  A-12.  Snapshot  Velocity  Vector  Plot 
at  Time  =  13.0  ysec 
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Figure  A-15.  Snapshot  Velocity  Vector  Plot  at  Time  »  6,66  psec 


Figu-r~  A-16.  Snapshot  Velocity  Vector  Plot  at  Time  -  7.19  nsec 


Figure  A-17.  Snapshot  Velocity  Vector  Plot  at  Time  =  7.78  Msec 


Snapshot  Velocity  Vector  Plot  at  Time  =8.71  psec 


Scale:  1  inch  (plotter  units)  =  u.t  cm 

1  inch  (plotter  units)  =  1.0  cm/ysec 
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Figure  A- 21.  Snapshot  Velocity  Vector  Plot  at  Time  =  9.68  psec 
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Figure  A- 24.  Snapshot  Grid  Plot  at  Time  =  5.95  ysec 


Snapshot  Grid  Plot  at  Time  =7.19  yisec 


Figuie  A- 28.  Snapshot  Grid  Plot  at  Time  =  8.31  psec 


Pigure  A-29.  Snapshot  Grid  Plot  at  Time  =  8.71  usee 
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Figure  A-30.  Snapshot  Grid  Plot  at  Time  =  °.14  ysec 


Figure  A-33.  Snapshot  Grid  Plot  at  Time  =  10.58  ysec 


APPENDIX  B 

GENERALIZED  BOUNDARY  POINT  MOTION 


88 


A  boundary  point  is  defined  as  any  non-interior  point  in 
the  grid  i.e.,  a  point  which  is  not  surrounded  by  four  zones. 

The  set  of  boundary  points  is  comprised  of  all  slavepoints, 
masterpoints ,  and  points  on  the  perimeter  of  the  grid. 

Possible  void  opening,  void  closing,  sliding  along  pistons, 
or.  masterlines,  and  the  various  combinations  of  these  options 
are  controlled  by  the  logic  shown  in  Figure  B-l  within  subroutine 
MOTION . 

B.l  SUBROUTINE  OPEN  (Figure  B-2) 

OPEN  performs  the  function  of  determining  if  a  point  has 
lifted  off  a  particular  segment.  The  calculation  is  a  simple 
one : 


(a)  Determine  the  distance  from  the  updated 
position  of  the  point,  P. 

(b)  Determine  if  this  distance  is  greater  than 
the  user-specified  distance  TOL.  If  it  is  not, 
the  point  does  not  lift  off;  return  to  subroutine 
MOTION.  If  it  is  greater,  the  point  may  possibly 
lift  off. 

(c)  Determine  for  the  case  where  the  distance 
is  greater  than  TOL  if  the  point  moved  into 
(crossed  over)  the  segment  or  moved  away  from 

the  segment.  If  the  point  moved  into  the  segment, 
the  void  does  not  open.  If  the  point  moved  away 
from  the  segment,  the  void  does  open. 
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Figure  B-l.  Logic  Controlling  Subroutine  NOTION 
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Figure  B-2  Subroutine  OPEN  * 

B.2  SUBROUTINE  SLIDE 

SLIDE  performs  the  function  of  moving  a  point  along  the 
geometric  constraint  of  a  piston  or  mas ter line  segment.  Friction 
may  or  may  not  be  present. 

The  flow  of  the  calculation  is  as  follows: 


a.  The  following  values  are  supplied  to  subroutine  SLIDE: 


XX,  YY  position  of  point  before  sliding 

XD,  YD  velocity  of  point  before  sliding 

XDD ,  YDD  acceleration  of  point  before  sliding 

b.  Locate  the  coordinates  and  velocities  for  the  line  segment 
(on  which  the  point  is  sliding) .  The  segment  may  be  defined  by 
two  points  in  space  (a  piston  segment)  or  by  two  masterline  points 
(a  master  segment) .  Slavepoints  may  slide  on  both  pistons  and 
masterline.  Other  points  may  only  slide  on  pistons. 
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Component  velocities  for  endpoints  1  and  2  are: 


XDl  = 

<XP1 

-  XZ1)/DLTH 

YD1  = 

(YP1 

-  YZD/DLTII 

XD2  = 

(XP2 

-  XZ2) /DLTH 

YD2  = 

(YP2 

-  YZ2) /DLTH 

where  DLTH  =  Atn+?s. 


c.  Resolve  the  point  acceleration  vector  into 
normal  and  tangential  components  with  respect  to 
the  segment  position  at  n  +  1,  PNDD,  and  PTDD. 

d.  Find  the  velocity  of  the  segment  at  point  P 
as  a  linear  weighted  average  of  the  velocities 
at  each  end  of  the  segment.  If  first  call  to 
SLIDE,  use  position  at  time  n  and  velocities  at 
n+J*.  If  not  the  first  call,  use  positions  at 

n  +  1  and  velocities  at  n+Jj.  This  difference 
between  the  first  and  subsequent  calls  is  to 
allow  for  the  point  to  slide  over  more  than  one 
segment . 


Segment  velocities  at  point  P(n+*s): 

SXDH  =  FRAC  *  XDl  +  (1.-  FRAC)  *  XD2 
SYDH  =  FRAC  *  YD1  +  (1.-  FRAC)  *  YD2 

Resolve  velocities  into  normal  and  tangential  components, 
and  STDH.  If  this  is  not  the  first  call,  set  SNDH  =  0.0. 


SNDH 
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e.  Resolve  point  P  velocity  into  normal  and 
tangential  components,  PND  and  PTD. 

f.  Compute  a  tentative  relative  sliding  velocity 
of  the  point  on  the  segment  at  n-fJj.  This  velocity 
is  for  use  in  calculating  a  possible  resisting 
tangential  force  due  to  friction.  Compute  resisting 
tangential  acceleration,  TDDRES. 

g.  Total  tangential  acceleration  of  point  P  is 
the  sura  of  its  free  and  resistive  components: 


TDD  =  PTDD  +  TDDRES 

h.  Update  normal  and  tangential  velocities  of  the 
point  n+i$ 


TDH  =  PTD  +  Atn  *  TDD 

NDH  =  SNDH 

Note  the  normal  velocity  of  the  point  is  identical  to  the  normal 
velocity  of  the  segment  at  point  P. 

i.  Find  new  position  of  the  point  P: 

XP  =  XX  +  Atn+3s  *  XD 

YP  =  YY  +  Atn+3*  *  YD 

j .  The  point  must  be  checked  to  see  if  it  still  lies 
on  the  segment  or  has  slid  off  one  of  the  ends.  The 
foilwing  possibilities  exist: 

(1)  The  point  remains  on  the  segment,  the 
motion  calculation  is  finished. 

(2)  Point  goes  off  an  end  of  the  segment 

»  If  there  is  not  another  segment 
connected  to  rhe  end,  the  point  may 

(a)  be  forced  to  stay  at  the  end  or 

(b)  be  allowed  to  become  a  free 
unconstrained  point.  The  user 
specifies  the  desired  option. 
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•  If  there  is  another  segment  connected 
to  the  end,  place  the  point  at  the  end 
of  the  segment  and  scale  the  accelera¬ 
tion  by  the  tangential  distance  moved. 
Return  control  to  subroutine  MOTION 
where  the  point  can  again  be  checked 
for  void  opening  if  desired.  If  no 
void  opens,  call  SLIDE  again  with  the 
resultant  acceleration  and  move  the 
point  along  the  next  segment. 


This  method  of  sequentially  moving  a  point  along  adjacent 
segments  defines  a  unique  and  accurate  updated  position  for  the 
point.  Previous  methods  of  utilizing  perpendicular  projections 
could  lead  to  non-unique  solutions  for  the  point  position  and 
possibly  large  anomalous  velocities. 

B.3  SUBROUTINE  CLOSE 

CLOSE  performs  the  function  of  determining  if  a  free  point 
has  impacted  a  particular  segment.  All  the  possible  segments 
are  checked  to  see  if  an  intersection  has  occurred.  For  slave- 
points,  both  masterline  segments  and  piston  segments  must  be 
checked.  Where  an  ambiguity  exists,  the  masterline  segments 
take  precedence  over  the  piston  segments.  For  a  regular  boundary 
point  including  a  masterpoint,  only  piston  segments  are  checked. 
The  loop  through  the  possible  segments  utilizes  the  following 
algorithm. 


A  quadrilateral  is  constructed  of  the  two  endpoints  of  the 
segment  at  times  n  and  n  +  1  in  a  frame  of  reference  where  the 
point  P  is  stationary. 
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An  automatic  rezone  package  has  been  developed  to  provide 
various  kinds  of  re zones  during  the  normal  processing  of  a 
Lagrange  2D  code .  The  unique  feature  of  the  system  is  that  it 
is  not  necessary  to  stop  the  problem  and  map  the  old  grid  onto 
the  new  grid  and  then  restart  the  calculation.  The  automatic 
rezoner  performs  its  rezoning  as  each  zone  is  calculated, 
according  to  the  options  desired. 

The  user  specifications  for  automatic  rezoning  are  as 
follows : 

(a)  Number  of  rezone  regions 

(b)  Indices  of  each  rezone  region,  IMNRZ  to  IMXRZ  and 

JMNRZ  to  JMXRZ 


For  every  rezone  region: 


(c)  NSTART 

(d)  NEND 

(e)  NFREQ 

(f)  N50PT 

(g)  NPO 


cycle  at  which  to  start  re zoning 

cycle  at  which  to  stop  rezoning 

frequency  of  cycles  at  which  re zones  are  to  be 
performed 

rezone  options 

printout  switch  for  re zone  check 
quantities 


Up  to  ten  re zone  regions  in  the  Lagrange  grid  can  be  speci¬ 
fied.  This  number  is  not  an  absolute  restriction  and  could  be 
easily  extended.  Various  rezone  regions  may  overlap  in  time  and 
space . 
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The  rezone  options  specified  in  the  packed  integer 
N50PT  =  N£N2N3N4N5  are  described  as  follows: 

N^  =  NHEZ ,  rezone  control 

NREZ  =  0  bypass  system 

NREZ  -  1  perform  rezone  of  zone  (JZ,IZ) 

NREZ  >  1  perform  NREZ  sweeps  (J  =  1, JMAJi )  on  column  IZ 


N2  =  NXY,  coordinate  re zoning 

NXY  =0  No  coordinate  re zoning 

NXY  =  1  Changes  the  old  coordinates  x0,  yQ  of  the 

point  being  re zoned  to  a  new  position  xp  and 
yp  where  xp  and  yp  are  constructed  by  a  mean 
parallelogram  method  using  neighbor  point 
positions. 

NXY  =  2  Performs  half  the  change  corresponding  to 
NXY  =  1.  That  is,  xnew  =  0.5 (Xn  +  xp)  and 
Ynew  =  0.5 (y0  +  Yp) •  This  option  is^ 
recommended  if  the  rezone  is  applied 
recurrently . 

NXY  =  3  Changes  the  old  coordinates  Xq,  yD  of  the 

point  being  rezoned  to  a  new  position  Xp  and 
yp  where  xp  and  yp  are  constructed  by  a 
centroid  method  using  neighbor  point 
positions . 

NXY  =  4  Performs  half  the  change  corresponding  to 
NXY  =  3.  That  is,  xnew  =  0.5 (Xq  +  Xp)  and 
Ynew  =  °*5(Yo  +  Yp) •  This  option  isF 
recommended  if  the  rezone  is  applied 
recurrently . 


N^  =  NMAP,  Interior  Map  Option 


The  rezoning  of  the  coordinates  of  a  point  will  result  in  a 
volume  change  of  up  to  four  zone  interiors.  Thus,  NXY  f  0 
should  be  accompanied  by  NMAP  f-  0.  The  interior  map  logic  con¬ 
serves  mass  density  and  energy  density  for  all  new  zones  whose 
volumes  have  been  reduced.  The  remaining  total  mass  and  inter¬ 
nal  energy  are  included  in  the  zones  that  are  expanded  according 
to  the  method  prescribed  by  the  value  of  NMAP,  If  the  total 
volume  of  the  rezoned  region  is  conserved,  then  NMAP  =  1  and  2 
give  identical  results.  Differences  result  when  a  material 
interface  or  free  surface  coordinate  re zone  has  resulted  in 
changed  volumes. 

NMAP  =  1  Conserves  mass  and  internal  energy;  stresses 
can  be  unrealistic. 

NMAP  =  2  Adds  or  subtracts  mass  and  internal  energy 
to  conserve  mean  stress. 


N4  =  NVEL,  velocity  rezoning 

NVEL  =  0-9  Indicates  the  weighting  factor  used  to 
weight  the  contribution  to  the  point's 
velocity  from  its  unrezoned  velocity  and 
the  velocities  of  its  neighbors.  The 
weighting  is  performed  according  to  the 
following  prescription:  neighbor 
velocities  are  weighted  by: 

NN 

1  ,  y  1  *  NVEL 

^ "  Ro  ti  -  RC  10 

where 

R  is  coordinate  position  of  point 
being  re zoned 

is  coordinate  position  of 
neighboring  points 

NN  is  the  number  of  neighboring 
points  up  to  8 
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The  original  velocity  is  weighted  by  the 
factor 


NVEL\ 

10  / 

Thus,  NVEL  =  0  applies  no  velocity  rezoning,  wher*  •*,  NVEL  =  9 
takes  90  percent  of  the  point's  velocity  from  surrounding  values 


Nj.  =  NMOM,  momentum  conservation 


Use  of  NMOM  >  0  will  provide  for  momentum  conservation  by 
remapping  of  the  velocity  of  the  re zoned  point  and  its  neighbors 


NMOM  =0  No  momentum  conservation. 

NMOM  =  1  Rigorous  conservation  of  mbmen.um 

independent  of  mass  conservation,  analogous 
to  NMAP  =  1.  Unrealistic  velocities  may 
result. 


NMOM  =  2  Conservation  of  momentum  density,  analogous 
to  NMAP  =  2. 


Certain  combinations  of  the  five  above  options  can 
obviously  lead  to  disaster.  An  intelligent  assessment  of  the 
objectives  of  a  calculation  should  be  made  when  applying  the 
automatic  rezoner  because  of  the  various  tradeoffs  of  a  particu¬ 
lar  option.  Reference  C.l  discusses  further  details  of  the 
automatic  rezoning  method. 


C.l  D.  E.  Maxwell,  An  Existing  Automatic  2D  Rezoner,  TCAM 

Technical  Note,  Physics  International  Company,  San  Leandro, 
California,  1973. 
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This  report  describes  a  technique  to  optimize  the  current 
shaped-charge  design  procedure  as  follows.  Starting  with  the 
desired  target  to  be  defeated,  a  determination  of  the  desired 
penetration  characteristics  of  the  jet  would  be  made.  Exist- 

trfL?e!ihPeneira?io"  thoory  would  then  be  U0ed  to  estimate  the 
ideal  characteristics  of  the  jet  to  defeat  the  given  taroet 

A  shaped  charge  launcher  would  then  be  designed^  give  these 
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ideal  jet  characteristics.  However,  a  suitable  design  pro¬ 
cedure  requires  (1)  a  viable  analytical  or  empirical  design 
approach  to  obtain  a  first  cut  shaped  charge  design,  (2)  a 
better  understanding  than  now  exists  of  the  detailed  mechan¬ 
isms  of  jet  formation,  and  (3)  a  better  understanding  of  the 
phenomenon  of  jet  penetration,  rhis  report,  which  is  con¬ 
tained  in  two  volumes,  addresses  the  first  two  of  these  re¬ 
quirements.  Volume  I  describes  the  use  of  the  existing  non- 
steady  state  theory  of  jet  formation  with  experimental  data 
and  one-dimensional  finite  difference  continuum  mechanics 
calculations  to  obtain  the  linear  collapse  velocity  for  gener¬ 
alized  axisymmetric  shaped  charges.  The  results  of  this  work 
are  then  used  to  obtain  non-unique  shaped  charge  designs 
which  give  the  required  idealized  jet  parameters.  Volume  II 
describes  the  modification  and  utilization  of  a  two-dimen¬ 
sional  finite  difference  continuum  mechanics  code  utilizing 
the  Lagrangian  coordinate  system  to  calculate  the  complete  jet 
formation  parameters  for  any  generalized  axisymmetric  shaped 
charge.  The  utilization  of  this  code  allows  a  more  detailed 
study  of  such  phenomena  as  jet  stability,  bifurcation  on  the 
axis,  shear  gradients,  viscosity,  shocks,  incipient  vaporiza¬ 
tion,  surface  tension,  and  possible  other  effects.  The  com¬ 
bined  use  of  both  the  engineering  formulations  along  with  the 
sophisticated  two-dimensional  code  calculation  allows  design 
engineers  the  versatility  to  design  the  most  optimum  shaped 
charge  for  their  particular  application. 
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